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Stochastic  inversion  is  a  well  known  technique  for  the  solution  of  inverse  problems  in  tomography. 

It  employs  the  idea  that  the  propagation  medium  may  be  represented  as  random  with  a  known 
spatial  covariance  function.  In  this  paper,  a  generalization  of  the  stochastic  inverse  for  acoustic 
travel-time  tomography  of  the  atmosphere  is  developed.  The  atmospheric  inhomogeneities  are 
considered  to  be  random,  not  only  in  space  but  also  in  time.  This  allows  one  to  incorporate 
tomographic  data  (travel  times)  obtained  at  different  times  to  estimate  the  state  of  the  propagation 
medium  at  any  given  time,  by  using  spatial-temporal  covariance  functions  of  atmospheric 
turbulence.  This  increases  the  amount  of  data  without  increasing  the  number  of  sources  and/or 
receivers.  A  numerical  simulation  for  two-dimensional  travel-time  acoustic  tomography  of  the 
atmosphere  is  performed  in  which  travel  times  between  sources  to  receivers  are  calculated,  given  the 
temperature  and  wind  velocity  fields.  These  travel  times  are  used  as  data  for  reconstructing  the 
original  fields  using  both  the  ordinary  stochastic  inversion  and  the  proposed  time-dependent 
stochastic  inversion  algorithms.  The  time-dependent  stochastic  inversion  produces  a  good  match  to 
the  specified  temperature  and  wind  velocity  fields,  with  average  errors  about  half  those  of  the 
ordinary  stochastic  inverse.  ©  2006  Acoustical  Society  of  America.  [DOI:  10.1121/1.2180535] 

PACS  number(s):  43.20.Dk,  43.28.We,  43.28.Vd  [AIT]  Pages:  2579-2588 


I.  INTRODUCTION 

Acoustic  tomography  is  widely  used  in  physics,  technol¬ 
ogy,  and  medicine  for  the  remote  sensing  of  inhomogeneous 
media.  When  applied  to  the  atmosphere,  acoustic  tomogra¬ 
phy  allows  one  to  estimate  (reconstruct)  temperature  and 
wind  velocity  within  a  tomographic  volume  or  area. 1-4  The 
practical  realization  of  acoustic  tomography  in  an  atmo¬ 
spheric  boundary  layer  was  reported  in  Refs.  5-10.  In  these 
tomography  experiments,  the  sources  and  receivers  of  sound 
were  located  on  masts  several  meters  above  the  ground  along 
the  perimeter  of  a  tomographic  area  that  was  a  square  or 
rectangular  with  side  lengths  of  several  hundred  meters.  The 
sound  travel  times  between  all  pairs  of  sources  and  receivers 
were  measured  and  used  as  input  data  for  inverse  algorithms 
to  estimate  temperature  and  wind  velocity  fields  within  a 
horizontal  slice.  This  kind  of  tomography  is  called  travel¬ 
time  tomography.  The  inverse  algorithms  used  were  the  sto¬ 
chastic  inversion  (SI)  approach5  and  the  simultaneous  itera¬ 
tive  reconstruction  technique. 6-10  In  the  present  paper,  a 
generalization  of  the  SI  approach  in  travel-time  acoustic  to¬ 


mography  of  the  atmosphere  is  developed  that  allows  one  to 
effectively  increase  the  number  of  data  without  increasing 
the  number  of  sources  and  receivers. 

The  idea  of  travel-time  acoustic  tomography  of  the  at¬ 
mosphere  is  based  on  the  fact  that  the  time  required  for 
sound  to  propagate  through  a  certain  volume  depends  on  the 
adiabatic  sound  speed  (and  hence  on  temperature)  and  wind 
velocity  within  that  volume.  More  specifically,  the  travel 
time  tj  of  the  sound  impulse  propagation  along  the  zth  ray 
can  be  expressed  as  the  following  integral  along  the  path  L, 
of  this  ray  (e.g.,  Ref.  11): 


Here,  Mg(R)  is  the  group  velocity  of  a  sound  impulse,  R  is 
the  position  vector  in  three  dimensions,  and  i=  1 ,2, ... ,/, 
where  I  is  the  number  of  ray  paths  of  sound  impulses  simul¬ 
taneously  propagating  through  this  volume  from  sources  to 
receivers.  It  can  be  shown1 1  that,  in  the  presence  of  wind,  the 
Ug  can  be  expressed  as  m^,  =  [cl(R)  + F2(R) 
+  2cL(R)n,(R)  •  V(R)]1/2.  Here  n,  is  the  unit  vector  normal  to 
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the  wave  front  of  a  sound  wave,  Y  is  the  vector  of  wind 
velocity,  and  cL  is  the  Laplace  adiabatic  sound  speed  that 
relates  to  the  acoustic  virtual  temperature  7av  by 

cl  =  yRaTm,  (2) 

where  y~  1.41  is  the  ratio  of  the  specific  heats  and  Ra  is 
the  universal  gas  constant  for  dry  air.  It  is  necessary  to 
take  into  account  the  specific  humidity  of  the  air  q  to 
obtain  the  values  of  the  actual  thermodynamic  temperature 
Tt  h  from  the  acoustic  virtual  temperature  by  using  the  fol¬ 
lowing  relationship:11  T.dv=  rth(l  +0.51  \q).  Note  that  the 
acoustic  virtual  temperature  differs  from  the  virtual  tem¬ 
perature  used  by  atmospheric  scientists.  When  deriving 
Eq.  (1),  it  is  assumed  that  ug  does  not  depend  on  time 
during  the  propagation  of  sound  impulses.  Since  in  the 
atmosphere  cL5>|V|,  this  assumption  is  valid  for  tomogra¬ 
phic  arrays  of  order  of  several  hundred  meters  considered 
in  Refs.  5-10  and  in  Sec.  IV.  The  goal  of  acoustic  travel¬ 
time  tomography  is  to  estimate  the  fields  Tav(R)  and  V(R) 
within  a  tomographic  volume,  given  the  travel  times  fj1 
and  the  locations  of  sources  and  receivers. 

There  are  different  techniques  available  to  solve  this 
problem.  Many  of  them  pursue  the  goal  of  finding  a  solution 
that,  after  substitution  back  into  Eq.  (1),  yields  travel  times 
as  close  to  the  measured  ones  as  possible.  These  methods  are 
quite  reliable  and  accurate  if  there  are  more  data  points  (e.g., 
travel  times)  than  unknown  model  values  (e.g.,  total  number 
of  points  at  which  the  fields  are  to  be  reconstructed).  How¬ 
ever,  this  condition  does  not  hold  for  atmospheric  tomogra¬ 
phy  unless  the  reconstruction  is  performed  with  low 
resolution.2'610’12,13  In  the  opposite  case,  when  the  number 
of  unknowns  is  greater  than  the  number  of  available  data 
points,  such  techniques  cannot  provide  a  unique  solution, 
and  some  additional  restrictions  must  be  imposed  on  the 
sought  fields.  Since  it  is  unknown  in  advance  whether  the 
actual  fields  being  estimated  obey  such  restrictions,  these 
techniques  may  yield  a  spurious  solution  that  perfectly 
matches  the  data  but  has  little  to  do  with  the  real  fields. 

In  contrast,  the  SI  approach  is  based  on  the  idea  that, 
using  the  available  data,  one  seeks  fields  that  have  the  mini¬ 
mum  average  deviation  from  the  real  ones.  This  approach 
also  requires  additional  information  about  the  sought  fields, 
namely,  the  sought  fields  are  treated  as  random  ones  with 
known  spatial  covariance  functions.1'3'5'14  Although  the  ac¬ 
tual  covariance  functions  in  the  turbulent  atmosphere  are  not 
exactly  known,  they  can  be  approximated  by  those  corre¬ 
sponding  to  the  von  Karman  or  Gaussian  spectra  of 
turbulence.111516  Note  that  other  techniques  that  partition  the 
fields  into  constant-valued  grid  cells  implicitly  assume  step¬ 
like  covariance  functions,  such  that  the  field  values  are  per¬ 
fectly  correlated  within  a  grid  cell  and  completely  uncorre¬ 
lated  outside.  Such  functions  are  much  less  plausible  than 
those  used  in  the  SI  approach.  The  applicability  and  advan¬ 
tages  of  the  SI  approach  in  travel-time  tomography  of  the 
atmosphere  with  high  resolution  are  demonstrated  in  Refs.  3 
and  5.  A  disadvantage  of  SI  is  that  it  does  not  use  the  fact 
that  the  sought  fields  are  correlated  not  only  in  space  but  also 
in  time. 


In  this  paper,  a  generalization  of  the  SI  approach  for 
acoustic  travel-time  tomography  of  the  atmosphere  is  devel¬ 
oped  that  allows  one  to  use  the  data  obtained  at  different 
times  to  reconstruct  the  temperature  and  wind  velocity  fields 
within  a  tomographic  volume  at  any  particular  time,  by  using 
spatial-temporal  covariance  functions  of  these  fields.  This 
generalization  is  called  time-dependent  stochastic  inversion 
(TDSI). 

The  general  idea  of  using  spatial-temporal  covariance 
functions  in  SI  is  known  in  the  literature.  For  example,  it  has 
been  successfully  used  in  satellite  altimetry  of  the  ocean  sur¬ 
face,  e.g..  Refs.  17  and  18.  In  altimetry,  the  deviations  of  the 
ocean  surface  height  from  its  average  level  are  known  at 
certain  times  and  spatial  points  (along  the  satellite  tracks). 
The  problem  is  to  estimate  these  deviations  at  other  points 
(in  space  and  time).  Thus,  this  problem  can  be  formulated  as 
the  space-time  interpolation  problem.  There  are  many  differ¬ 
ent  techniques  to  interpolate  the  data.  However,  SI  interpo¬ 
lates  them  in  such  a  way  that  the  covariance  of  the  resulting 
fields  holds  true.  The  application  of  SI  for  these  purposes 
was  proposed  in  Ref.  19. 

In  contrast  to  satellite  altimetry,  travel-time  tomography 
cannot  be  formulated  as  an  interpolation  problem,  since  the 
fields  that  are  subject  to  estimation  are  unknown  everywhere. 
Therefore,  the  mathematical  formalism  of  TDSI  in  acoustic 
travel-time  tomography  of  the  atmosphere  differs  from  that 
in  altimetry.  In  the  present  paper,  this  mathematical  formal¬ 
ism  is  developed.  To  verify  that  the  TDSI  approach  improves 
the  quality  of  the  reconstruction,  a  numerical  experiment  was 
carried  out  for  two-dimensional  (2-D)  travel-time  acoustic 
tomography  of  a  horizontal  atmospheric  layer. 

The  paper  is  organized  as  follows.  In  Sec.  II,  a  general 
theory  of  2-D  TDSI  is  developed  and  a  formula  for  the  op¬ 
timal  stochastic  inversion  matrix  is  derived.  Calculations  of 
the  covariance  matrices  that  appear  in  the  inversion  matrix 
are  presented  in  Sec.  III.  In  Sec.  IV,  the  numerical  experi¬ 
ment  of  2-D  travel-time  acoustic  tomography  of  the  atmo¬ 
sphere  is  described.  Some  aspects  of  the  SI  and  TDSI  algo¬ 
rithms  are  discussed  in  Sec.  V.  The  conclusions  are  presented 
in  Sec.  VI. 


II.  THEORY 

In  this  section,  the  theory  of  2-D  TDSI  is  developed. 
There  are  three  steps  that  precede  this  development.  First,  it 
is  taken  into  account  that  Tm  and  V  depend  not  only  on  the 
spatial  coordinates  but  also  on  time  t.  Therefore,  tj  depends 
on  time  t  as  well.  Second,  Eq.  (1)  is  linearized  due  to  specific 
conditions  in  the  atmosphere,  and,  with  the  same  degree  of 
accuracy,  the  travel  paths  L ,  are  approximated  by  straight 
lines.  Third,  it  is  shown  how  to  reconstruct  the  spatial  mean 
values  of  the  temperature  and  wind  velocity  fields  within  a 
tomographic  area  with  the  help  of  the  conventional  least 
square  estimation.  Finally,  the  problem  for  the  TDSI  ap¬ 
proach  in  the  reconstruction  of  temperature  and  wind  veloc¬ 
ity  fluctuations  is  stated  and  a  general  solution  of  this  prob¬ 
lem  is  obtained. 


2580  J.  Acoust.  Soc.  Am.,  Vol.  119,  No.  5,  May  2006 


Vecherin  eta/.:  Time-dependent  stochastic  inversion 


A.  Linearization 

Let  u( r . t)  and  u(r,t)  be  x  and  y  components  of  the 
two-dimensional  vector  of  wind  velocity  V: 

V(r,f)  =  u(r,t)ex  +  v(r,t)ey.  (3) 

Here  t  is  time,  a  2-D  vector  r  specifies  a  spatial  point  within 
a  tomographic  area  with  the  Cartesian  coordinates  (x,y),  and 
ec  and  ev  are  the  unit  vectors  along  the  x  and  y  axes,  respec¬ 
tively.  The  adiabatic  sound  speed  cL,  temperature  7’.lv,  and 
wind  velocity  fields  within  a  tomographic  area  at  time  t 
can  be  represented  as  sums  of  their  spatial  average  values 
Co(r),  T0{t),  uo(t),  and  u0M  and  their  fluctuations  c(r,f), 
T(r,t),  u(r,t),  and  v(r,t ): 

cL(r,t)  =  c0(t)  +  c(r,t),  Tjr,t)  =  T0(t )  +  T{r,t), 

(4) 

u(r,t)  =  u0(t )  +  u(r,t),  u(r,f)  =  u0(t)  +  v(r,t) . 

Since  in  the  atmosphere  the  absolute  values  of  the  adiabatic 
sound  speed  fluctuations  and  the  wind  components  are  much 
smaller  than  the  spatial  average  value  of  c0,  Eq.  (1)  can  be 
linearized  to  the  first  order  of  these  fluctuations  and  the 
travel  paths  of  sound  can  be  approximated  by  straight 

r  5,11 

lines: 

,tr/ a  Li  (  ,  Mo(f)cos  (fj  +  u0(f)sin 

^ - J 

+  u(r,t)sin  (p,j ,  (5) 

where  1=1,2, ... ,/  denotes  the  ray’s  path  number,  7  is  the 
total  number  of  travel  paths,  L,  is  the  length  of  the  ith  ray 
path,  r  £  Lh  and  tp,  is  the  angle  of  the  ith  ray  relative  to  the 
positive  direction  of  the  x  axis.  Note  that  the  noise  s,  in 
the  1th  measured  travel  time  was  added  to  the  right-hand 
side  of  Eq.  (5).  Furthermore,  when  deriving  this  equation 
the  following  relationship  between  the  fluctuations  of 
temperature  and  adiabatic  sound  speed  was  used:  c(r,f) 
=  [co(r)/27’0(r)]7’(r,t).  This  relationship  follows  from  Eqs. 
(4)  and  (2),  with  accuracy  up  to  the  first  order  of  these 
fluctuations.  The  inverse  problem,  which  is  studied  in  this 
paper,  is  to  reconstruct  cq,  T0,  T,  n,  and  v,  given  tf,  Lh  and 
Vi- 

B.  Reconstruction  of  mean  fields 

The  mean  values  T0,  u0,  and  v(}  within  a  tomographic 
area  at  any  time  can  be  reconstructed  with  the  help  of  the 
least  square  estimation  by  using  the  travel  times  obtained  at 
the  same  moment  of  time.  For  this  purpose  one  should  set 
the  fluctuating  parts  of  these  fields  to  zero,  so  that  the  inte¬ 
gral  in  Eq.  (5)  vanishes.  Then  Eq.  (5)  can  be  rewritten  in 
matrix  notation: 

Gf  =  b,  (6) 

where  the  elements  of  the  column  vector  b  are  known  and 
given  by  i>,  =  r'r(f) /L„  the  unknown  column  vector  f  has 


three  elements,  /,  =  l/c0(r),  f2=u0(t)/cl(t),  f3  =  v0{t)  I  c20{t), 
and  the  matrix  G  is  given  by 


G  = 


1  -  cos  <p j  -  sin  <p[ 


1  -  cos  i pj  -  sin  i pi 


(7) 


Then,  using  the  7  elements  of  the  vector  b  (7  is  assumed 
to  be  greater  than  3),  one  should  solve  the  overdetermined 
problem  for  the  three  unknowns  (the  elements  of  the  vector 
f)  using  the  least  square  estimation: 

f^G^'G^b,  (8) 


extract  from  them  the  values  of  c0,  n0,  and  u0,  and  calculate 
T0  using  Eq.  (2). 


C.  Time-dependent  stochastic  inversion 

Once  the  spatial  mean  values  of  temperature,  adiabatic 
sound  speed,  and  wind  velocity  fields  within  a  tomographic 
area  are  known,  it  is  worthwhile  to  introduce  the  column 
vector  of  data  d(t)  obtained  at  time  t  with  elements 

dj(t)  =  L,(c0(t)  -  u0{t) cos  <p;  -  u0(f)sin  <P;]  -  cl(t)tf(t)  +  £(f) , 

(9) 

where  noise  includes  the  errors  of  travel  time  measure¬ 
ments  Sj(t)  and  errors  in  the  estimation  of  70,  c0,  u0,  and  u0. 
Using  Eq.  (5),  this  data  vector  can  be  expressed  as  d(t) 
=d0(t)  +  ^(t),  where  the  elements  of  the  noise-free  data  col¬ 
umn  vector  da(t)  are  given  by 

d0iU)  =  J  d/^^^T(r,t)  +  u(r,t)cos<pi  +  v(r,t)sin(p^. 

(10) 

The  problem  now  is  to  reconstruct  the  fields  of  fluctua¬ 
tions  ^(r^o),  «(r,f0),  and  u(r,70)  at  any  chosen  time  f0. 
Thus,  the  column  vector  of  models  at  this  time  moment  t0  is 
given  by 

m(f0)  =  [7’(r1,r0);  ...  ;T(r y,t0); 


n(r[,f0);  ...  ;«(ry,f0) ;u(r!,r0) ;  ...  ;u(ry,f0)],  (11) 

where  J  is  the  number  of  spatial  points  within  the  tomogra¬ 
phic  area  at  which  the  fields  are  being  reconstructed;  here 
and  in  what  follows,  the  semicolon  between  elements  de¬ 
notes  that  these  elements  are  arranged  in  a  column. 

Let  us  assume  that  one  can  collect  data  (fV+1)  times 
during  the  time  Nt  by  performing  sound  scans  of  the  tomog¬ 
raphic  area  and  forming  the  vector  d  at  each  scan.  Suppose 
that  the  last  scan  was  performed  at  the  time  t,  which  can  be 
earlier,  equal  to,  or  later  than  the  time  t0  at  which  one  tries  to 
reconstruct  the  models.  In  this  paper  it  is  assumed  that  these 
scans  are  performed  at  equal  time  intervals  r,  although  such 
an  assumption  is  not  necessary.  The  temperature  and  wind 
velocity  fields  in  this  tomographic  area  are  changing  with 
time,  which  yields  different  data  vectors  at  each  scan.  In  this 
case,  at  the  time  t  there  will  be  data  available  from  (77+1) 
scans.  That  is,  at  the  moment  t  one  can  form  the  vector  d  of 
all  available  data: 
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(12) 


(14) 


d  =  [d(t  -  Nr)-,d(t  -  Nt+  t);  ...  ;<7(r)]. 


Here  d  is  a  column  vector  containing  the  ( /V  + 1 )  column- 
vectors  d(t-nr),  each  of  length  7,  where  n=Q,l, ...  ,N. 
Therefore,  the  total  length  of  the  data  vector  d  is  (IV +  I )7. 


1.  Problem  statement 

Using  the  data  vector  d  given  by  Eq.  (12),  it  is  necessary 
to  find  a  linear  estimation  m(t0)  of  unknown  models  m(r0)  at 
time  t0.  The  case  f0  >  1  corresponds  to  advanced  estimation; 
the  case  t0<t  corresponds  to  posterior  estimation;  and  t0=t 
corresponds  to  a  real  time  reconstruction.  In  all  these  cases, 
one  can  use  available  data  obtained  from  all  scans  to  estimate 
the  unknown  fields.  If  one  uses  only  the  data  obtained  at 
moment  t0  to  estimate  the  models  m(r0)  at  the  same  moment, 
that  would  correspond  to  standard  SI  in  atmospheric  travel¬ 
time  tomography.5 

In  TDSI,  the  estimation  of  the  unknown  fields  is  sought 
in  the  same  form  as  in  standard  SI: 

m(r0)  =  Ad,  (13) 

where  unknown  elements  ajk  of  matrix  A  must  be  deter¬ 
mined.  Here  j=  1 ,2, ...  ,37,  where  3/  is  the  length  of  mod¬ 
els,  and  k=  1 ,2, . . .  ,(N+ 1)1.  We  now  introduce  the  column 
vector  of  discrepancy  e  between  the  true  and  reconstructed 
fields  at  time  t0: 


R,nd  =  [Vmd(to,t-NT),BnJt(),t-NT+  r),  ...  ,Bmrf(f0,r)], 


R, 


dd  ' 


B dd(t  ~  Nr,  t  -  Nr)  B dd(t  -  Nr,  t  -  Nt  +  t) 

B dd(r~  Nr+  r,t  —  Nt)  B dd(t-Nr  +  T,t  —  Nr+  t) 


Kddit’t-Nr) 


B  dd{t,t-Nr+  r) 


=  ’hjito)  -  ntj{t0) . 

To  find  the  elements  a/7f,  let  us  require  that  they  give  the 
minimum  of  the  elements  of  the  mean  square  errors  vector 

(e2): 

(dj) l->  min(e?),  (15) 

iajk} 

where  parentheses  (■)  denote  the  mathematical  expectation. 
These  mean  square  errors  are  the  diagonal  elements  of  the 
error  covariance  matrix  liee: 

R«-(«7),  (16) 

where  e  is  the  column  vector  of  discrepancy  given  by  (14). 

2.  Solution 

The  matrix  A  that  solves  the  problem  stated  above  is 
given  by  the  same  formula  as  for  ordinary  SI,  since  the  deri¬ 
vation  of  this  formula  does  not  depend  on  a  particular  struc¬ 
ture  of  the  models  m  and  data  d  (see  Appendix  A): 

A  =  RmdRdd>  (17) 

where  Rmd  =  (mdr)  and  Rdd  =  (dd7)  are  model-data  and 
data  covariance  matrices. 

Since  the  data  and  models  have  been  chosen,  as  shown 
in  Eqs.  (12)  and  (11),  these  matrices  have  the  following 
block  structure: 


Kdd(t-Nr,t) 
®dd(t-Nr+  r,t) 


(18) 


(19) 


Here  Blm/(f , ,  t2)  =  (m(ti)dT(t2))  is  the  covariance  matrix 
of  size  [3/,/]  between  the  models  at  time  moment  t1  and 
data  at  moment  t2,  Bdd(t1,t2)  =  (d(tl)dT(t2))  is  the  covariance 
matrix  of  size  [/,/]  between  data  at  moment  tl  and  data  at 
moment  t2,  and  the  comma  between  the  elements  of  the  Rmd 
matrix  denotes  that  the  Bmd  matrices  are  arranged  in  a  row. 
Thus,  the  dimensions  of  the  Rmd  and  Rdd  matrices  are 
[3/,(A/+l)7]  and  [(N+  1)1, (N+ 1)7],  correspondingly.  Note 
that  the  Rdd  matrix  is  symmetric. 

It  is  assumed  that  noise  in  the  data  £(f)  is  independent  of 
the  sought  fields  and  the  noise-free  data  vectors,  is  not  cor¬ 
related  with  itself  for  different  time  moments  and  for  differ¬ 
ent  paths,  and  can  be  described  by  the  normal  distribution 
N( 0,O|).  Therefore,  in  the  presence  of  noise. 


Rmd-Rmd0>  (20) 
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Rdd  -  Rd0d0  +  °^I,  (21) 

where  I  is  the  identity  matrix.  Since  noise  in  the  data  can 
easily  be  taken  into  account  by  these  formulas,  further  con¬ 
siderations  will  be  focused  on  the  noise-free  data  d0. 

It  can  be  shown  (see  Appendix  B)  that  the  error  covari¬ 
ance  matrix  Re6  that  corresponds  to  such  an  estimation  of 
models  [see  Eq.  (13)]  with  the  optimal  matrix  A  given  by  Eq. 
(17)  is 

R«re  =  Rmm  -  RmdRddl^md  >  (22) 

where  Rmm  =  (in(?(l)m7( r0))  is  the  models  covariance  matrix 
at  time  f0.  The  elements  of  the  main  diagonal  of  the  Ree 
matrix  are  equal  to  expected  mean  square  errors  of  the  re¬ 
construction  at  each  point.  They  can  be  calculated  without 
knowledge  of  the  original  fields.  These  averaged  errors  can 
differ  from  the  errors  of  the  reconstruction  of  a  particular 
realization  of  the  original  fields. 
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Note  that  the  developed  mathematical  formalism  of 
TDSI  is  valid  for  statistically  nonstationary  and  inhomoge¬ 
neous  random  fields. 


III.  COVARIANCE  MATRICES 

Since  the  optimal  stochastic  inverse  operator  A  given  by 
Eq.  (17)  is  determined  in  terms  of  the  Rmd  and  Rdd  matrices, 
it  is  worthwhile  to  consider  some  important  particular  cases 


when  these  matrices  can  be  calculated  explicitly.  These  con¬ 
siderations  will  be  carried  out  taking  2-D  travel-time  tomog¬ 
raphy  of  the  atmosphere  as  an  example. 

A.  2-D  travel-time  tomography 

In  the  case  of  travel-time  tomography,  one  can  use  the 
linear  relationship  between  d0  and  the  sought  fields  given  by 
Eq.  (10)  to  obtain  an  expression  for  the  covariance  matrix 
Bmd(|(f| ,  t2)  between  the  models  at  time  tl  and  the  noise-free 
data  at  time  t2: 


[Bmdo(fi,f2)]ji  =  {mj(tt)d0t(t2)) 


■l 


co(h) 

dl\  r,°  ,]  Smjih)T{rd2))  +  (mj{ti)u(r,t2))cos  <pf  +  <w/f1)i;(r,t2))sin  <p; 


2 T0(t2) 
Coih) 


JL. d'{  2T0(t2) f  1 ; r’ ?2)  j  ’  if  1 

dl(Buu(rj,tl;r,t2)cos(pi  +  Buv(rj,tl;r,t2)sm(pi),  ifj+l^;^27,  (23) 

dl(Bvu(rj,t1;r,t2)cos  (pt  + B^r^; r,f2)sin  <p,),  if  27  +  1  «  37, 

where  i=l,2,...,7,  /=  1 ,2, ...  ,37,  r  eL(,  Btt,  Buu,  Bvv,  Buu,  and  Bvu  are  the  spatial-temporal  covariance  functions  of  the 
corresponding  fields  marked  as  the  subscripts,  and  the  ry  are  the  chosen  spatial  points  within  the  tomographic  area  at  which  the 
sought  fields  are  reconstructed;  these  points  stay  fixed  during  the  integration. 

Similarly,  an  expression  for  the  covariance  matrix  Bd  d  {t\  ,t2)  between  the  noise-free  data  at  time  r,  and  time  t2  is  given 
by 

=  (dOi(ti)doP(h))=  \  dl\  dl'\  /'^fl^!2\grr(r7i;r',f2)  +  fi,,„(r,ti;r',f2)cos  y,cos  ipp 


+  5„„(r,f1;r',t2)sin  cpi  sin  +Buu(r,t1;r',f2)cos  sin  ip  +B„„(r,r1;r',f2)cos  <p  sin  <p,-  f , 


(24) 


where  i,p=l,2, ... ,/,  r  e  L„  r'  e  Lp.  Note  that 
Bvu(r , ft ; r' , f2) -Buv(r’  ,t2-,r,tl).  When  deriving  Eqs.  (23) 
and  (24),  it  is  assumed  that  BIu  =  BIv  =  0. 

Once  the  matrices  Bmrf()(f  x ,  t2)  and  B d^Jj\ ,  h)  are  calcu¬ 
lated,  one  should  use  Eqs.  (18)  and  (19)  to  form  the  Rmd0  and 
Bd0d0  matrices  and,  then,  Eq.  (21)  to  take  noise  into  account. 


B.  Stationary  fields 

If  the  temperature  and  wind  velocity  fields  are  statisti¬ 
cally  stationary,  their  covariance  functions  depend  only  on 
the  difference  of  their  temporal  arguments: 

B(r,ti  ;r',f2)  =  B(r.r’,At),  A  t=t2  —  t\.  (25) 

The  notation  B  without  subscripts  stands  for  the  covariance 
function  of  any  two  fields  considered  above.  In  this  case, 
Eqs.  (18)  and  (19)  are  modified  as  follows: 
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Rmd  =  [BmrfO  -Nt- t0)  ,Bmd(f  -Nt+  T-t0),  ...  , 


Bmrf(f-fo)L  (26) 

Brfd(O)  B  dd(r)  ■■■  B  dd(Nr) 

D  B  dd(-r)  Bdd(0)  •••  B  m(Nt-t) 

Kdd  -  •  •  •  • 

_B  dd(-NT)  B  dd(-NT+r)  Bdd(0) 

(27) 

Note  that,  for  any  At, 

Brfd(Ar)  =  B>Tdd(-  Ar).  (28) 


It  is  sufficient  to  know  only  the  first  cell  row  of  the  matrix 
Rdd  since  it  has  the  same  matrices  along  each  diagonal  (a 
block  Toeplitz-like  structure).  Also  note  that  the  Rdd  matrix 
is  independent  of  time  r0,  i.e.,  for  the  given  set  of  data  it  is 
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necessary  to  calculate  this  matrix  only  once  to  reconstruct 
fields  at  any  time.  Another  interesting  fact  is  that  the  matri¬ 
ces  Brfrf(0)  on  its  main  diagonal  are  equal  to  the  covariance 
matrix  Rdd  of  ordinary  SI. 

C.  Frozen  turbulence 

If  turbulence  may  be  considered  frozen,  each  spatial 
point  of  the  temperature  and  wind  velocity  fields  is  moving 
with  a  constant  speed  U.  In  this  case,  the  temperature  field  at 

time  t2  can  be  expressed  in  terms  of  the  field  at  time  tl  by  the 

20 

following  relationship: 

T(r,t2)  =  T(r-VAt,tl)-  (29) 

The  same  relationship  is  valid  for  the  turbulent  wind  velocity 
fields.  This  allows  one  to  obtain  the  spatial-temporal  covari¬ 
ance  matrices,  which  are  unknown  in  general,  in  terms  of  the 
spatial  covariance  matrices.  The  latter  can  be  modeled  by  the 
covariance  functions  corresponding,  for  example,  to  the  von 
Karman  or  Gaussian  spatial  spectra  of  turbulence.  For  ex¬ 
ample,  for  the  temperature  field,  we  have 

BTT{r,tx-y  ,t2)  =  (r(r,ri)r(r',f2)) 

=  (7’(r,f1)T(r'-U  At,tl))  =  Bsri(r,r'  -U  At), 

(30) 

where  BSTT  is  the  spatial  covariance  matrix  of  the  temperature 
field.  The  last  equality  is  true  for  statistically  stationary 
fields.  Similar  equations  can  be  derived  for  the  spatial- 
temporal  covariance  functions  of  the  other  fields.  Since  fro¬ 
zen  turbulence  is  a  particular  case  of  stationary  fields,  Eqs. 
(26)-(28),  stay  valid  for  this  case  as  well.  Thus,  for  frozen 
turbulence  the  spatial-temporal  covariance  matrices  Rmd  and 
Rdd  can  be  expressed  in  terms  of  known  or  inferred  spatial 
covariance  matrices. 

The  frozen  turbulence  hypothesis  is  widely  used  in  me¬ 
teorology  and  the  study  of  wave  propagation  in  random  me¬ 
dia,  e.g.,  Refs.  20  and  21.  This  hypothesis  is  valid  with  a 
good  accuracy  for  eddies  in  the  inertial  range  of  turbulence. 
However,  the  applicability  of  this  hypothesis  to  eddies  in  the 
energy  range  is  questionable.  Note  that  in  the  present  pa¬ 
per,  the  frozen  turbulence  hypothesis  is  used  to  concretize 
the  spatial-temporal  covariance  functions.  TDSI  can  still  be 
used  if  this  hypothesis  is  invalid. 

IV.  NUMERICAL  EXPERIMENT 

In  this  section,  a  numerical  experiment  that  demon¬ 
strates  the  significant  improvement  in  the  quality  of  tempera¬ 
ture  and  wind  velocity  field  reconstructions  possible  with 
TDSI  is  described.  A  snapshot  of  the  temperature  and  wind 
velocity  fields  was  created  by  large  eddy  simulation  (LES),23 
which  is  widely  used  in  meteorology  (e.g..  Refs.  24  and  25). 
LES  numerically  solves  the  Navier-Stokes  equations  for  ed¬ 
dies  large  enough  to  be  resolved  by  the  spatial  grid.  The 
effect  of  smaller  eddies  on  the  resolved  flow  is  parametrized. 
LES  produces  realistic  wind  velocity  and  temperature  fields 
for  flow  processes  that  are  well  resolved  by  the  grid. 

Using  the  single  snapshot  of  the  LES  temperature  and 
wind  velocity  fields,  the  perfectly  frozen  turbulent  fields  with 


x(m) 


FIG.  1 .  (Color  online)  The  layout  of  sources  and  receivers  in  the  numerical 
experiment. 

the  velocity  components  I/t=4  m/s  and  f/v=0m/s  were 
created.  The  travel  times  were  calculated  each  second  (t 
=  1  s)  using  Eq.  (5).  Two-dimensional  linear  interpolation 
was  employed  to  obtain  the  value  of  the  integrand  in  Eq.  (5) 
along  the  travel  paths;  then  the  integral  was  calculated  nu¬ 
merically.  These  travel  times  were  disturbed  by  white  noise 
with  standard  deviation  crE=5  /xs,  which  corresponds  to  rela¬ 
tive  errors  of  order  1%  to  42%  in  the  data  d0(t).  The  total 
number  of  such  numerical  scans  was  18,  and  the  problem 
was  to  estimate  the  temperature  and  wind  velocity  fields 
rav(r,f0),  w(r,f0 ),  and  u(r ,r0)  at  time  t0=9  s.  The  time  origin 
coincided  with  the  first  scan.  The  number  of  sources  as  well 
as  receivers  was  5,  so  that  the  number  of  rays  I  at  each  scan 
was  25.  The  tomographic  area  was  a  square  of  80  by 
80  meters.  The  layout  of  acoustic  sources  and  receivers  and 
the  corresponding  ray  paths  are  shown  in  Fig.  1. 

The  reconstruction  consisted  of  two  stages.  In  the  first 
stage  the  spatial  mean  values  of  the  temperature  and  wind 
velocity  fields  T0(f0),  m0(/0),  and  u0(f0 )  within  the  tomogra¬ 
phic  area  were  reconstructed  using  the  approach  described  in 
Sec.  II  B.  The  true  and  estimated  spatial  mean  values  of  tem¬ 
perature  and  wind  fields  are  presented  in  Table  I. 

One  can  see  from  Table  I  that  the  reconstruction  of  spa¬ 
tial  mean  fields  was  very  accurate,  the  differences  between 
true  and  estimated  mean  values  was  0.14  K  for  the  tempera¬ 
ture  field,  0.03  m/s  for  the  x  component  of  the  wind  velocity 
vector,  and  0.01  m/s  for  its  y  component. 

In  the  second  stage,  the  fluctuations  of  the  temperature 
and  wind  fields  7(r,f0),  u(r,t0),  and  v(r,t0)  from  their  spa¬ 
tial  mean  values  were  reconstructed  using  either  standard  SI 
(only  data  at  time  t0  were  used,  the  total  number  of  data 


TABLE  I.  Actual  and  estimated  spatial  mean  values  of  temperature  and 
wind  velocity  fields. 


Mean  fields 

T0  (K) 

«0  (m/s) 

Vo  (m/s) 

True 

301.73 

3.09 

1.73 

Estimated 

301.87 

3.06 

1.71 

2584  J.  Acoust.  Soc.  Am.,  Vol.  119,  No.  5,  May  2006 


Vecherin  et  al.\  Time-dependent  stochastic  inversion 


x(m) 


x  (m) 


x  (m) 


FIG.  2.  The  original  and  reconstructed 
temperature  (K)  and  wind  velocity 
(m/s)  fields  of  fluctuations  at  time  mo¬ 
ment  r0=9  s.  (a)  Original  u  field,  (b) 
Original  v  field,  (c)  Original  T  field, 
(d)  The  u  field  reconstructed  by  SI.  (e) 
The  v  field  reconstructed  by  SI.  (f) 
The  T  field  reconstructed  by  SI.  (g) 
The  u  field  reconstructed  by  TDSI.  (h) 
The  v  field  reconstructed  by  TDSI.  (i) 
The  T  field  reconstructed  by  TDSI. 


points  was  equal  to  the  number  of  rays  /=  25)  or  TDSI  (data 
obtained  at  18  different  times  were  used,  which  resulted  in 
450  data  points).  To  describe  the  spatial  covariance  of  the 
temperature  and  wind  velocity  fields  within  the  tomographic 
area,  the  following  Gaussian  covariance  functions  were 
used:11 


r,r')  =  oy  exp 


(r-r')2 


(31) 


(r,r')  =  (Tu 


exp 


r\2 


(r-r')2)/,  (y-y') 


,  (32) 


®uu(r>r')  =  (Tv  exp 


,,  2  /  (r-r')-)/,  (x-x') 


'\2 


(33) 


„s  /  r\  I  (r  -  r')2\  (x- x')(y -y') 

®uy(r,r  )  =  (rHav  exp)  - - y -  - y - ,  (34) 


l- 


l 


where  crT,c rlncrv  are  the  standard  deviations  of  the  corre¬ 
sponding  fields,  lT  and  /  are  their  correlation  lengths,  r 
=  (x,y),  and  r'  =  (x',y').  Then,  the  relations  similar  to  that 
given  by  Eq.  (30)  were  used  to  obtain  the  spatial-temporal 
covariance  functions  of  the  sought  fields.  The  noise-free  ma¬ 
trices  Rmdo  and  Rd0d0  were  formed  with  the  help  of  Eqs.  (26) 
and  (27),  where  the  Bmd0  and  B matrices  were  calculated 
using  Eqs.  (23)  and  (24);  noise  was  taken  into  account  by  Eq. 
(21). 

In  the  described  numerical  experiment,  there  were  five 
parameters  in  the  SI  and  TDSI  algorithms  that  must  be  cho¬ 


sen:  oy,  cr,(,  <x„,  lT.  and  /.  These  parameters  were  estimated 
from  the  original  LES  fields:  crr=0.14  K,  cr„  =  0.72  m/s,  au 
=  0.42  m/s,  and  lT=l=  15  m. 

The  original  and  reconstructed  fields  were  described  by 
matrices  of  size  21  by  21.  After  all  calculations  were  done, 
the  fields  were  interpolated  (2-D  linear  interpolation)  for  il¬ 
lustrative  purposes.  The  fluctuating  parts  of  the  original 
fields  are  presented  in  Figs.  2(a)-2(c).  Figures  2(d)-2(f) 
show  the  reconstructed  fields  using  the  SI  approach.  The 
reconstruction  results  of  TDSI  are  given  in  Figs.  2(g)-2(i). 
As  one  can  see,  standard  SI  matched  the  contours  of  the 
sought  u  field  fairly  well,  but  the  v  and  T  fields  were  recon¬ 
structed  poorly.  In  contrast,  the  time-dependent  stochastic 
approach  allowed  the  much  more  detailed  and  accurate  re¬ 
construction  of  all  fields.  To  characterize  the  expected  im¬ 
provement  of  the  reconstruction  it  is  worthwhile  to  consider 
the  expected  errors  of  the  reconstruction  given  by  Eq.  (22).  It 
is  convenient  to  normalize  these  mean  square  errors  by  the 
corresponding  field  variances  so  that  their  values  lie  in  [0,  1], 
If  these  normalized  mean-squared  errors  (NMSE)  are  of  or¬ 
der  unity,  the  errors  of  reconstruction  are  of  the  order  of 
variance  of  the  original  field,  i.e.,  one  has  a  poor  reconstruc¬ 
tion.  Conversely,  if  NMSE  are  zeros,  it  is  expected  that  the 
reconstructed  and  actual  fields  are  identical.  NMSE  can  be 
recalculated  into  root  mean-square  errors  (RMSE).  NMSE 
for  SI  are  presented  in  Figs.  3(a)-3(c)  while  Figs.  3(d)-3(f) 
show  NMSE  for  TDSI.  To  characterize  the  overall  quality  of 
the  reconstruction,  NMSE  at  each  spatial  point  and  corre¬ 
sponding  RMSE  were  averaged  over  the  tomographic  area. 
These  averaged  values  of  NMSE  and  RMSE  are  presented  in 
Table  II. 
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FIG.  3.  (Color  online)  The  expected 
normalized  mean  square  errors,  (a) 
The  u  field  reconstructed  by  SI.  (b) 

The  v  field  reconstructed  by  SI.  (c) 

The  T  field  reconstructed  by  SI.  (d) 

The  u  field  reconstructed  by  TDSI.  (e) 
The  v  field  reconstructed  by  TDSI.  (f) 
The  T  field  reconstructed  by  TDSI. 


It  follows  from  Fig.  3  and  Table  II  that,  for  TDSI,  the 
average  NMSE  is  about  four  times  less  for  the  T  and  v  fields 
and  seven  times  less  for  the  u  field  than  the  corresponding 
NMSE  for  ordinary  SI.  This  corresponds  to  a  46%  reduction 
of  the  RMSE  for  the  T  field,  a  52%  reduction  for  the  v  fields, 
and  a  62%  reduction  for  the  u  field. 

V.  DISCUSSION 

The  idea  of  the  stochastic  approach  is  based  on  the  as¬ 
sumption  that  covariance  matrices  of  the  sought  fields  are 
known.  That  means  that  in  the  case  of  the  Gaussian  covari¬ 
ance  functions  all  five  parameters  (variances  and  correlation 
lengths  of  the  sought  fields)  should  be  known  in  advance.  In 
the  practical  implementation  of  travel-time  tomography, 
however,  these  parameters  are  unknown.  One  can  estimate 
these  parameters  by  measuring  the  temperature  and  wind  ve¬ 
locity  fluctuations  with  the  use  of  conventional  meteorologi¬ 
cal  sensors  or  using  a  turbulence  similarity  theory. 

Furthermore,  the  covariance  of  the  isotropic  turbulence 
in  the  atmosphere  is  better  described  by  the  von  Karman 
spectrum  rather  than  by  the  Gaussian  spectrum.  The  latter 
was  used  in  the  numerical  experiment  because  it  was  then 
possible  to  calculate  the  Rllld  matrix  analytically  and  use 
only  a  single  numeric  integration  to  get  the  Rdd  matrix.  This 


TABLE  II.  Average  expected  errors  of  the  reconstruction  of  temperature  and 
wind  velocity  fields  by  SI  and  TDSI. 


Fields 

T 

u 

V 

Errors 

NMSE 

RMSE  (K) 

NMSE 

RMSE  (m/s) 

NMSE 

RMSE  (m/s) 

SI 

0.96 

0.13 

0.56 

0.53 

0.63 

0.33 

TDSI 

0.24 

0.07 

0.08 

0.20 

0.15 

0.16 

significantly  improved  the  accuracy  and  speed  of  the  calcu¬ 
lations.  Since  the  actual  fields  were  not  described  by  the 
Gaussian  covariance  functions,  the  reconstruction  quality  im¬ 
provement  seems  rather  unexpected.  However,  there  is  a 
relatively  simple  explanation  of  this  fact.  In  a  certain  spatial 
scale  range,  the  covariance  functions  of  atmospheric  turbu¬ 
lence  can  be  approximated  by  the  Gaussian  covariance  func- 

27 

tions  with  appropriate  variances  and  correlation  lengths. 
Therefore,  the  Gaussian  covariance  functions  can  be  used  for 
SI  and  TDSI  as  approximations  of  the  actual  ones.  It  is  ex¬ 
pected  that  the  use  of  the  covariance  function  corresponding 
to  the  von  Karman  spectrum  of  turbulence  will  improve  the 
reconstruction. 

One  of  the  advantages  of  the  stochastic  approach  is  that 
noise  in  the  data  plays,  in  a  certain  sense,  a  positive  role  by 
regularizing  the  Rdd  matrix  [see  Eq.  (21)].  Indeed,  the  pres¬ 
ence  of  noise  just  adds  an  additional  term  to  the  main  diag¬ 
onal  of  the  Rdd  matrix,  which  improves  its  condition  (keeps 
it  invertible),  although  it  does  smooth  the  solution  somewhat. 
This  regularization  is  especially  important  for  TDSI,  where 
the  condition  of  the  matrix  Rdd  may  be  very  poor  and  may 
lead  to  spurious  solutions. 

Figures  2(a)-2(i)  show  that  the  reconstruction  of  the 
temperature  fluctuation  field  was  not  as  good  as  that  for  the 
velocity  fluctuations.  The  probable  reason  for  this  is  a  small 
effect  of  the  temperature  fluctuations  on  the  travel  times  in 
comparison  to  the  wind-velocity  fluctuations.  Indeed,  the 
LES  temperature  fluctuations  were  in  the  range 
[-0.15,0.35]  K  while,  for  many  meteorological  problems, 
the  acceptable  error  of  temperature  measurements  is 
±0.3  K.  ’  Note  that  to  improve  the  reconstruction  of  tem¬ 
perature  fluctuations  in  acoustic  tomography,  one  can  use 
reciprocal  transmission  of  sound  waves.8'26 
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The  proposed  TDSI  algorithm  also  allows  one  to  esti¬ 
mate  the  temporal  mean  values  of  the  temperature  and  wind 
velocity  fluctuations  at  any  given  spatial  point.  This  might  be 
important  when  the  temporal  changes  of  these  fluctuations 
are  negligible  during  the  time  interval  Nt.  For  this  purpose,  it 
is  necessary  to  set  all  temporal  arguments  in  the  matrices 
Rm(]  and  Rdd  equal  to  zero.  Note  that  in  this  case  the  noise 
variance  cr|  should  be  enlarged,  because  the  noise  includes 
not  only  the  measurement  uncertainty  but  also  the  neglected 
variation  of  the  original  fields  during  the  observed  time  in¬ 
terval. 

VI.  CONCLUSION 

In  this  paper,  a  generalization  of  the  stochastic  inversion 
approach  for  tomographic  problems  in  the  atmosphere  was 
developed.  The  key  idea  of  this  generalization  is  that,  by 
using  spatial-temporal  covariance  functions  of  temperature 
and  wind  velocity  fields,  one  can  use  data  obtained  at  differ¬ 
ent  times  to  reconstruct  these  fields  at  any  particular  time. 
This  allows  one  to  significantly  enlarge  the  effective  amount 
of  tomographic  data,  while  keeping  the  total  number  of 
sources  and/or  receivers  fixed.  The  efficiency  of  the  devel¬ 
oped  method  was  demonstrated  in  a  numerical  experiment 
that  simulated  the  acoustic  2-D  travel-time  tomography  prob¬ 
lem  of  the  atmosphere.  This  numerical  experiment  showed  a 
remarkable  improvement  of  the  reconstruction  when  one 
used  TDSI  instead  of  standard  SI. 
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To  find  the  elements  ajk  of  matrix  A,  let  us  require  that 
they  give  the  minimum  of  each  element  of  the  mean  square 
errors  vector  (e2): 

(^)^mm(e2),  (A3) 

{«/*} 

where  k=  1,2,...  ,D. 

To  satisfy  the  requirement  (A3),  one  should  take  the 
partial  derivative  along  the  current  element  ajp  of  the  matrix 
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“ — ((ajkdk  -  nij)2)  =  2((alkdk  -  m^d  ) 
daip 

=  2 (aik(dkdp)  -  (midp)) 

=  2(a//r.[Rdd]jt/,  -  [Rmd]i;,)  =  0,  (A4) 

where  the  equality  ilajkl dalp  =  SjjSpk  was  taken  into  account, 
and  the  S’ s  are  Kronecker’s  delta  symbols.  The  last  equation, 
after  canceling  the  factor  2,  can  be  written  in  matrix  notation: 

ARdd  -  Rmd  =  0 .  (A5) 

Therefore,  under  assumption  that  the  inverse  matrix  Rdd  ex¬ 
ists,  the  optimal  matrix  A  is  given  by 

A  =  RmdRdd  ■  (A6) 


APPENDIX  B:  ERROR  COVARIANCE  MATRIX 
DERIVATION 


This  material  is  partly  based  upon  work  that  was  sup¬ 
ported  by  the  U.S.  Army  Research  Office  under  Contracts 
No.  D AAD 19-03-1-0104  and  No.  DAAD19-03-1-0341. 
Also,  we  would  like  to  thank  E.  G.  Patton  and  P.  P.  Sullivan 
of  the  National  Center  for  Atmospheric  Research  for  provid¬ 
ing  the  LES  data.  The  LES  was  supported  in  part  by  the  DoD 
High-Performance  Computing  Modernization  Program.  Fi¬ 
nally,  we  would  like  to  thank  anonymous  reviewers,  whose 
comments  allowed  us  to  improve  this  paper. 

APPENDIX  A:  OPTIMAL  STOCHASTIC  OPERATOR 
DERIVATION 

Let  m  and  d  be  two  arbitrary  random  column-vectors  of 
lengths  M  and  D,  correspondingly,  with  a  known  matrix  of 
the  second  statistical  moment  Rmd  =  (md7).  Furthermore,  let 
the  matrix  Rdd  =  (dd7)  also  be  known.  Note  that,  in  the  case 
when  (m)  =  0,  (d)  =  0,  the  Rmd  and  Rdd  are  the  covariance 
matrices.  The  problem  is  to  estimate  the  unknown  vector  m 
if  a  particular  realization  of  the  random  vector  d  is  given.  To 
do  this,  one  can  seek  for  the  estimation  m  of  the  unknown 
vector  m  in  the  form 

rii  =  Ad.  (Al) 


If  a  random  column  vector  r/  relates  to  another  random 
column  vector  £  by  a  known  nonrandom  matrix  C  as  r] 
=  C  f ,  then  its  covariance  matrix  R  nv=(  w')  is  expressed  in 
terms  of  the  covariance  matrix  R^=  (^r)  as 

R^=CRffCr  (Bl) 


In  the  case  of  travel-time  tomography,  the  column  vector  of 
discrepancy  is  given  by  e=  rh-m,  where  m  =  Ad.  Therefore, 
the  discrepancy  can  be  expressed  in  terms  of  the  matrix  mul¬ 
tiplication  using  block  matrix  notation: 


6=  [A -I] 


d 

m 


(B2) 


where  I  is  the  identical  matrix.  Using  the  formula  (Bl),  one 
can  get  the  expression  for  covariance  matrix  R££: 


R« 


[d7!!!7] 


(dd7)  (dm7) 


At 

-I 


(md7)  (mm^JL-I. 


-  ARddAr-  AR7d  -  RmdA7  +  Rimn.  (B3) 


Introduce  the  column  vector  of  discrepancy  e  between 
the  true  and  reconstructed  vectors: 

€j  =  riij-  m.j,  (A2) 

where  j=  1 ,2, ...  ,M. 


Note  that  this  formula  for  the  covariance  matrix  of  the  dis¬ 
crepancy  vector  is  valid  for  any  matrix  A.  Since  the  optimal 
matrix  A  of  the  stochastic  inversion  and  time-dependent  sto¬ 
chastic  inversion  is  given  by  A  =  RmdRdd,  the  R££  matrix 
takes  the  form 
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Ree  -  I^rnm  _  ^md^dd^md  ■  (B4) 

Note  that  the  diagonal  elements  of  the  error  covariance  ma¬ 
trix  Reer  represent  the  mean  square  errors  of  the  estimation 
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